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A SQUARE ROOT METHOD FOR COMPUTATIONALLY EFFICIENT MODEL 

PREDICTIVE CONTROL 

[0001] This invention was conceived in performance of work under U.S. 
Government Contract N00421-01-2-0131. 

BACKGROUND OF THE INVENTION 

[0002] This invention relates to generally to control systems and more 
particularly to model predictive control systems. 

[0003] Model predictive control logic has many advantages in controlling 
practical multivariable systems. A multivariable system has multiple goals (also known 
as commands) and has multiple effectors that change system dynamics, for achieving 
those goals. Multivariable systems may also include significant cross-coupling, in which 
key effectors each significantly effect multiple goals. Therefore, a controller for this 
system should also be multivariable to decouple the responses. In a decoupled system, 
when a goal changes, the control algorithm responds by setting multiple effectors so that 
this one goal changes without significantly changing the other goals. Furthermore, the 
system cross-coupling is often dynamic, changing in time, where some coupling is 
persistent, some transitory, some fast, some slow, or it may first couple in one direction, 
then reverse and couple in the other. Therefore the control algorithm decoupling should 
also be dynamic, evolving in time to counter the changing nature of the coupling. In 
addition to de-coupling, the multivariable control logic must also provide the tradition 
feedback properties of reducing the effects of changes in the system dynamics and the 
effects of disturbances, all while ensuring the system remains stable, safe, and durable. 
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[0004] There is much control theory, some reduced to practice, supporting the 
design of feedback control logic for coupled multivariable systems. Most of this theory 
presumes the system to be controlled has no inequality constraints, referred to as limits 
here, on any of the system variables or effectors. This is a very significant mismatch 
between theory and practice. Real and practical systems have substantial and significant 
limits. Some limits are physical, such as effector ranges, and some are operational, to 
ensure stability, safety, and durability. In a real sense, if no limits are hit, the system has 
been over designed. This mismatch has led to great difficulties and ad hoc methods in 
control law design, including splitting logic into many modes, where there is a separate 
control mode to cover each permutation of limit or fault condition, and more ad hoc logic 
to dynamically control the switching between separately designed dynamic modes. 
Coordinating these modes and ensuring there are no adverse dynamic interactions 
between modes (e.g., repeatedly jumping back and forth between two modes, sudden 
jumps, etc.) requires a significant amount of additional ad hoc logic. Ad hoc design 
makes the control software difficult to design, difficult to maintain and greatly reduces 
the degree to which the software can be reused on a similar system. 

[0005] Model predictive control (MPC) is a multivariable control theory that 
explicitly includes limits and, therefore, provides a good match with practical systems. 
MPC can also be configured to respond in realtime to changes in the system, such as 
actuator faults. Thus, MPC provides a formal method of control algorithm design for 
multivariable systems that can decouple responses, as well as physically possible, even as 
limits are hit and faults occur. In MPC design, there is no need for ad hoc designs of 
extra control modes to handle limit conditions and faults. This is expected to provide 
significant increases in controller software design productivity, make it easier to upgrade 
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and otherwise maintain control software, and make it possible for the control system to 
accomplish more difficult tasks. The latter include more autonomous operation, more 
completely integrating system responses, and operating with reduced physical stresses. 

[0006] The most significant roadblock to implementing model predictive 
control is the amount of computation required. That is why the initial practice of MPC, 
which began 20 to 30 years ago in the chemical industry, was for systems with relatively 
slow dynamics. A chemical plant controller updates actuator settings on the order of 
once every 5 minutes. These early applications stressed the computer power available at 
the time. In a gas turbine engine, the update rate must be around 40 times per second, an 
approximately 12,000-fold difference, and the control function may be more complicated. 
Even so, MPC may be practical for engines, and other vehicle applications, because 
computers are now available that are 500 to 1000 times faster. That still leaves a 10 to 
30-fold gap, depending on system complexity, between the computation MPC appears to 
require and computer power available. This invention is critical because it provides 
numerical techniques for reducing the computational burden of model predictive control 
by 20 to 100-fold, closing the gap. 

[0007] Model Predictive Control 

[0008] Model predictive control theory has been reduced to practice, most 
notably in the chemical industry. In support of this, model predictive control theory has 
been described in numerous technical papers and textbooks. It is summarized here. The 
main objective of this sunmiary is to show that implementing MPC requires a large 
amount of onboard computation and that approximately 90% of this computation 
involves numerous solutions of a large matrix equation in realtime. This section shows 
that when the optimality equations and variables are ordered in a certain way, the matrix 
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equation has a particular sparse, banded, symmetric structure. These properties imply the 
matrix includes many zeros, the nonzero entries are gathered near the diagonal, and a 
matrix "square root" exists. This particular form is developed in this section and termed 
the Large Sparse Matrix Equation (LSME) equation. The square root method exploits 
this structure in a novel way to greatly reduce computation. 

[0009] In model predictive control, actuator conmiands, u, are determined as 
the result of solving an optimization problem in realtime. Model predictive control is a 
type of digital controller. Digital controllers are based on a computer that periodically 
samples operator commands and measured responses. As a consequence of each sample, 
the digital controller specifies new actuator settings, which is referred to as a controller 
update. These settings are most often held until the next update. MPC specifies these 
settings as the first time point of an actuator trajectory. The trajectory specifies u(n), for 
controller update times numbered with integer n, where n varies from L to N+L-1. That 
is, for N time points. The current sample time is enumerated L. This actuator trajectory 
is determined as the unique trajectory that minimizes the performance index, PI, where 



[0010] and where the summation is over the N future time points, x^ is a 
vector of the states of a system dynamic model and u„ are the actuator commands at each 
of the N time points. The matrices C and D are coefficients of linearized output 
equations, which relate system outputs for which there are goals to the dynamic model 
states and actuator settings. Q and R are lower triangular weighting matrices, and vectors 
f and g are driving terms that are functions of a desired trajectory for the system output, 
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desired or nominal actuator values, and weighting matrices Q and R. Lower triangular 
matrix M is a terminal-weighting matrix and may be present for a variety of control 
theoretic reasons. 

[0011] The nominal system and, therefore, the optimization must obey the 
state dynamics model of the system 

[0012] which optimization procedures term equality constraints. The 
optimization must also obey the physical and operational limits linked to the 
Limited variables yc(n) = Cc^ • + Dc„ 'U^-\-a^ <Y^ 

[0013] for the N time points. The first equality represents a model of how the 
limited variables are linked to the state dynamics model. The explicit inclusion of the 
limit equation in the theory and resulting logic design is a key advantage of MPC. 

[0014] A MPC designer specifies the control algorithm by determining the 
state dynamics model, the limit model, and the performance index. Performance index 
design involves selecting the variables to include in output vector, y, and choosing the 
weighting matrices Q, R, and M. The control software designer must also provide 
reusable software for solving the optimization problems in realtime. 

[0015] The above optimization problem is a quadratic programming problem: 
a special type of constrained mathematical optimization. Quadratic programming is 
much more reliable than general purpose constrained optimization, and the associated 
restrictions are mild. The two major methods for solving quadratic programming (QP) 
problems on a computer are Active Set (AS) methods and Interior Point (IP) methods. 
The invention described below will significantly reduce the amount computation needed 
with either. 
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[0016] In MPC operation, the state vector at the current sample time, , is 
given. It is a function of measurements of the system response up to the controller update 
time L. Thus, the MPC controller is a feedback controller. The net output of the above 
optimization is the vector u^, the first time point in the u trajectory, which is the actuator 
command for the L* update of the MPC controller. 

[0017] The same MPC optimization must be performed, with the variables 
shifted in time, to determine the next actuator conmiand, u^^^. The MPC optimization 
must be re-solved every controller update period. The next three sections summarize 
how the above type of optimization is most conmionly solved in an analytical fashion. 

[0018] Ad joint methods for Absorbing Constraints Into J 

[0019] In the usual method of solution for the optimal control, the equality 
constraints (the state dynamics model) are adjoined to the performance index with 
Lagrange multipliers, pn, and the inequality constraints (the model of system variables 
with respect to their physical limits) are adjoined with multipliers, mn, to form the 
augmented performance index as follows. 
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and the terminal weight Pn = ]i^N • 



[0020] In the Lagrange method, the adjoined constraint terms are formulated 
so that the constraints are satisfied when the augmented performance index is optimized 
as if there were no constraints. The adjoined terms also have zero value at the optimal. 
Thus adjoining the constraints does not change the value of optimal J. For inequality 
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constraints the adjoint method also requires that for each limit, at each time n, either the 
corresponding m=0 or the inequahty constraint is on its limit. 
[0021] Karush-Kuhn-Tucker (KKT) conditions 

[0022] A set of algebraic equations whose solution specify the optimal system 
trajectory are known as the (KKT) conditions. For the above MFC optimization, it is 
known that the solution includes the following KKT conditions for each controller update 
between time L and time N+L-1: 
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Cc„-x„+Dc,.«„+a„+f„-y, =0 

the split boundary conditions : 
Xj^=xfeedback and Pn+l = M^+l^a^+l^n+l 

[0023] and the additional conditions that all t's and m's must be greater than 
zero. If these equations are solved, the actuator commands, u(n), that are part of the 
solution, used as inputs to state equation model of the system dynamics, will minimize 
the performance index, while meeting the limits. 

[0024] General Large Sparse Matrix Equation for Quadratic 
Programming 

[0025] The KKT conditions for the N time points can be collected into the 
following large sparse matrix equation (LSME). This must be repeatedly solved when 
solving the MFC quadratic programming problem, using either Interior Foint (IF) or 
Active Set (AS) methods. The LSME has the form: 
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[0026] where z is a vector containing the state, x, control, u, and state 
equation adjoint variable, p, for each time point, group by time, m is the adjoint variables 
for the inequality constraints, grouped by time, f and K are vectors. H and J are banded 
matrices, where 
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[0027] T is a diagonal matrix and, the matrices are shown for the controller 
update at time L=0. Matrix H is also symmetric, which implies there is a square root 
matrix, Hr, such that Hr'*Hr = H. The LSME has the same form at other time points but 
with different parameter values. The LSME is common to both of the major quadratic 
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program solvers, except that the value of T and the right hand side (RHS) vector depend 
on the quadratic program solver used. With either method the LSME must be solved 
repeatedly (10 to several hundred times, depending on system complexity) for each 
controller update. 

[0028] This form of the LSME is novel and special. It is achieved by 
ordering the optimality equations (rows of the LSME) and the variables in z in a special 
way. Then the resulting LSME matrix is sparse, banded, and symmetric. 

[0029] Interior Point Methods 

[0030] In the Interior Point method for solving quadratic programs, T is 
diagonal. Each element is the square root of the product of the corresponding slack 
variable, t, and adjoint variable, m, for each time point. J is sparse banded matrix, and 
includes all the limit equations. For the most popular family of interior point methods, T 
and the right hand side vector vary twice each iteration. The right hand side vector is 
determined from the solution of the previous iterate. H and J remain constant for the 
iterations at the same controller update period. The solution for the optimal actuator 
settings will require 10 to several hundred iterations, depending on problem complexity. 
For more explanation, see Ref. 1, especially their case where the slack variable has been 
eliminated. 

[0031] Active Set Methods 

[0032] The Active Set method searches for which limits are active. Only 
some of the inequality constraints, Jc*z <= Ymax, will be active. Each iteration of an 
Active Set qpsolver "guesses" which of the constraints are active. Let the vector "active" 
contain the indices of the active constraints. When active, the inequalities become 
equalities, and Jc(active,:)*z = Ymax(active), and the system is running on the active 
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limits. Compared to the LSME, J = Jc(active,:), K = Ymax(active), and T=0. In a given 
controller update period, a number of iterations of "active" will be evaluated in a search 
for the correct active set. The search strategy determines a new iteration of "active" 
based on the results of the previous iterate' s solution of the LSME. H will remain 
constant for the iterations for a given controller update, 

[0033] The key observation here is that some form of the LSME must be 
solved numerous times to solve the MPC quadratic program each controller update, 
whether interior point or active set solver is used. For practical controllers, the LSME is 
solved 10 to 100 times per controller update, depending on N and other factors. 



SUMMARY OF THE INVENTION 

[0034] The novel algorithm described here improves upon this prior art by 
providing a square root method for factoring the Large Sparse Matrix Equation. Solving 
matrix equations such as the LSME involves first factoring the matrix into the product of 
an upper triangular matrix, one or more diagonal matrices, and a lower triangular matrix 
(or vice versa) and then solving easier matrix equations, two of which involve triangular 
matrices and the rest diagonal matrices. 

[0035] Among other novel features described in more detail below, the 
present invention uses a novel method to form LSMroot directly, without the need to ever 
form the Large Sparse Matrix itself. The present invention also uses a novel method to 
form Hr directly once per controller update, without the need to first forrri H, and uses Hr 
throughout the quadratic program solution. 
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[0036] These and other novel features of the present invention reduce the 
computational burden of model predictive control by 10 to 30-fold. This reduction in the 
computational burden makes the use of model predictive control practical for more 
complicated systems with much faster dynamics, such as a gas turbine engines. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[0037] Other advantages of the present invention can be understood by 

reference to the following detailed description when considered in connection with the 

accompanying drawings wherein: 

[0038] FIGURE 1 is a schematic of a model predictive control system 

according to the present invention. 

[0039] FIGURE 2 is a flowchart of the method according to the present 

invention utilized by the system of FIGURE 1. 

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT 
[0040] The algorithm of the present invention provides a square root method 
for factoring the Large Sparse Matrix Equation. Solving matrix equations such as the 
LSME involves first factoring the matrix into the product of an upper triangular matrix, 
one or more diagonal matrices, and a lower triangular matrix (or vice versa) and then 
solving easier matrix equations, two of which involve triangular matrices and the rest 
diagonal matrices. For the expected MFC problem sizes, the factorization is one to two 
orders of magnitude more complex than the solution of the diagonal and both triangular 
matrix equations combined. This invention is an algorithm for performing this 
factorization much faster. 
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[0041] At the highest level of description, this algorithm achieves great 
numerical efficiency by: 

[0042] 1. Ordering the equations and variables to get a synmietric, banded 

LSME. 

[0043] 2. Realizing this LSME needs to be solved many times each 
controller update 

[0044] 3. Realizing that because the LSME matrix, LSM, is symmetric that it 
has a complex square root matrix, LSMc such that LSM= LSMc'*LSMc and LSMc is 
lower triangular. 

[0045] 4. Generalizes this concept of a square root matrix to use only real 
numbers (for computational efficiency) so that LSM = LSMroot'*S* LSMroot, S is a 
constant diagonal with only 1 or -1 entries along the diagonal, and LSMroot is lower 
triangular. 

[0046] 5. Uses a novel method to form LSMroot directly, without the need to 
ever form LSM. 

[0047] 6. Realizing that at a given controller update, the matrix H remains 
constant throughout the quadratic program solution and that it has a square root, Hr. 

[0048] 7. Using a novel method to form Hr directly, without the need to 
form H. This involves a sequence of smaller well-known QR factorizations. 

[0049] 8. Using a novel method to form LSMroot for each iteration of the 
quadratic program solver by solving one of two generalized QR factorizations, each of 
which reuse Hr. 

[0050] 9. Novel algorithms for performing the generalized QR factorization. 



12 



Express Mail No. EV221421154US Docket No. 67097-006gEH- 10871 

[0051] 10. Using standard methods for completing the solution of the LSME 
once the LSM has been factored into form of Step 4. 

[0052] These steps embody two key novel means for the efficient 
factorization of the LSME. The first is the efficient factorization of the submatrix H, 
which only needs to be done once per controller update (Steps 1, 6, and 7). The second is 
reuse of the H factorization to more quickly factor the entire LSME (Steps 8 and 9), 
which needs to be done 10 to 200 times per controller update (once or twice per QP 
iteration), depending on the QP method and size of problem. 

[0053] Once per Controller Update: Factor H 

[0054] Input the submatrices of H and J: A. B. C, D, Q*'^, R''', and Cc and 
Dc, which describe the MPC problem for this controller update. Factor H as follows. 

[0055] Step I: Compute the N separate and smaller QR factorizations, using a 
standard QR factorization algorithm 
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[0056] Start with n=N and work backwards in time to N-l, N-2, 1. The 
left hand side matrix is the input to the factorization; the right hand side is the result. In 
the standard methods, it is optional whether the orthogonal matrix, Ortho, is produced. It 
is not needed here. The matrices A, B, C, D, Q^'^ and R"^ are specified ahead time, as is 
M^. The first QR factorization then determines Mj^.j, which is an input to the second 
factorization for time N-l, and so on. Thus all N can be performed recursively, backward 
in time. This recursive sequence of QR factorizations and their value in solving the 
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unconstrained MPC problem (if J and K equal 0) are known in control theory textbooks. 
Their use in factoring H (Step 2) is novel. As is the reuse of factored H in manner that 
allows H to be factored only once per controller update (Steps 3, 4, and 5). 

[0057] Step 2: Form matrices L and W where L is lower triangular and W is 
block diagonal, as follows. Also, form and store the diagonal elements of matrix S, 
whose entries are 1 for except for columns corresponding the adjoint variables, p, in 
vector z, in which case the entries are -1. 



'0 

0 / -M 



L = 



-I 



-I 



-I 
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' N-l 

0 / 



-M 



(2) 



(3) 



[0058] Standard methods for efficient matrix computation would store the 
same information in L, S, and W in a more compact manner. The sparse matrix form is 
used here to aid comprehension. 
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[0059] There is no need to determine H in this algorithm. However, if H is 
desired it can be computed with 

[0060] H = Hr'*S*Hr, where Hr = L W\ (4) 

[0061] To save computations, Hr is not explicitly formed either. Instead, L 
S, and W are separately used. 

[0062] Every Quadratic Program Solver Iteration: Factor and Solve a 

LSME 

[0063] Input the parameters that describe the next LSME that needs to be 
solved from the QP solver search method: T, K, and f. In general, quadratic program 
solution involves a structured, iterative search for the solution. The search method and 
the parameters used to specify the LSME are different with AS and IP. Accordingly, the 
remaining steps of the square root algorithm differ for IP and AS methods. 

[0064] Step 3(IP): Perform the following generalized QR factorization, using 
the generalized QR algorithm described in the next section. 



= U ' Pin, where Pin = 



(5) 



[0065] The matrix in square brackets on the right is the input; the other 
matrices are the output. U does not need to be computed here. P is a lower triangular 
matrix. Since T is diagonal and J and L are banded, P will be banded as well and the 
special QR factorization can be used to exploit the banded structure to reduce 
computations. This is a novel, generalized factorization because U in not orthonormal. 
Instead, U has the more general property 
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[0066] The matrix (termed LSM) in the LSME does not need to be form here. 
However, if desired, it can be computed as follows. 
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[0067] Step 4(IP): Solve the intermediate equation 
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[0068] for wl and w2 via a block back substitution comprising the substeps: 

4. 1 . Solve wl = K, for w2, using a standard method for diagonal matrices 

4.2. Solver •ql = w2 forql, using a standard method for diagonal matrices 

4.3 Solve q2 = W • f - J 1 , using a method that exploits the structue of W and J 

4.3. w 1 = ql for wl , using a standard method for backsubstitution for banded 
matrices 

[0069] Step 5(IP): Solve the final equation for mz and m. 



(8) 



[0070] via the block forward substitution method using the substeps: 

5.1,Solveq3 = S • wl, using a method that exploits the structureof S 
5.2 Solve P m2 = 93 for mz, using a standard method for forward substitution 
for banded matrices 

5.3. Solve q4 = T"^ J • mz, for q4, using a method that exploits the structure of T and J 

5.3. Solve T • m = w2 - q4 , for m, using a standard method for diagonal matrices 

5.4. If the variables, p, are needed : solve z = W^mz, using a standard method for 

back substitution, modified to exploit the structure of W (the p's are the only elements of z that differ from mz^ 
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[0071] The LSME is now solved. The variables m and z (or mz) are output 
to the search method of the IP method. 

[0072] Step 3(AS): Perform the following generalized QR factorization, using 
the generalized QR algorithm described in the next section. 

5y'P = 5*''t/ Pm, where Pin = [L''^y^J (9) 

[0073] The matrix in square brackets on the right is the input; the other 
matrices are the outputs. U does not need to be computed here. P is upper triangular. 
Unfortunately, the banded structure of is lost when it is inverted, thus, the QR 
factorization cannot exploit the banded structure to save computation and P is not banded. 
This is a novel, generalized factorization because U is not orthonormal. Instead, U has 
the more general property 

U'^S'U = S 

[0074] L V is an upper triangular matrix, except for a band of nonzero entries 
below the diagonal. The generalized QR factorization should exploit the fact that the 
lower triangular part below this band is already zero. 

[0075] There is no need to determine the matrix LSM. However, it can be 
computed as follows. 
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[0076] 
[0077] 



Matrix T is zero in the AS method. 

Step 4(AS): Solve the intermediate equation 
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[0078] for wl and w2 via block back substitution using the substeps: 

4. 1 . Solve iJwl^W'f, for w2, using a standard method for back substitution 

4.2. Solve ^1 = 5 • w2 , using a method that exploits the structure of S 

4.2. Solve L- ^2 = ^1 for q2, using a standard method for forward substitution 

4.3. ^3 = (K + J • q2), using a method that exploits the banded structure of J 

4.4. wl -q?> for 1 , using a standard method for forward substitution 

[0079] Step 5(AS): Solve the final equation 



[0080] for m and mz via the block forward substitution method using the 
substeps: 

5.1. Solve qr4 = -5^ • w 1 , using a method that exploits the structure of Sr 

5.2 Solve P 'm = q4 for m, using a standard method for back substitution 

5.3 Solveq5 = /^m, using a method that exploits the structure of J 

5.4 Solve L^q6 = <75, for q6, using a standard method for back substitution 
6. Solve ql = 5(w2 - q6) , using a method that exploits the structure of S 

6. Solve L • mz = ^7 for mz, using a standard method for forward substitution 

7. If the variables, p, are needed : solve z = W^mz, using a standard method for 
back substitution, modified to exploit the structure of W 

[0081] Generalized QR Factorization Algorithm 

[0082] A novel generalized QR factorization is required in Step 3 (AS or IP). 
Well known algorithms are available to determine a traditional QR factorization of the 
form P = U*Pin, where Pin is input matrix, U has the property that U^U=I, and P is either 
upper triangular or lower triangular. The methods described are a novel extension of a 
family of known QR methods using a sequence of Householder transformations. 

[0083] The Step 3 factorizations in the square root method require 
Sr^^^P = S^^^U * Pin , where Pin and S are inputs, Sr and S are diagonal matrices whose 




(12) 
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diagonal entries are either 1 or -1, P is lower triangular or upper triangular, Sr is an 
output, and U has the property that U^S*U=S. Note that if Pin is real, then P is also real. 
U may contain complex number components. This factorization can be performed in the 
same manner as the well-known method of QR factorization using a sequence of 
Householder transformations, if the following transforms are used in place of the 
Householder transforms. 



[0084] and e^ is a column vector of zeros except for a 1 entry in the ith 
location. Xj and e^ are defined following the same procedure as in the traditional QR 
factorization. The above description using matrices and vectors is for explanatory 
purposes only. Those skilled in computation linear algebra know that the algorithm that 
implements this could save computations by not actually forming U and v. Those so 
skilled also know how to arrange the sequence of transformations, define each x based on 
the input matrix, and apply the transformations to produce either a lower triangular or 
upper triangular output matrix, P. The output matrix, Sr, is a diagonal matrix where each 
element is the product of the corresponding sw and S elements. As a further explanatory 
aid, when S is equal to the identity matrix, all sw's will be one and the new transform is 
reduced to the traditional Householder transform. Other than the substitution for the 
Householder transform, tradition QR factorization algorithms can be used to produce the 
factorizations in Step 3 (AS or IP) and can be used to exploit the structure, which is 
optional, in Step 3 IP. Note that in Step 3 IP, the structure of Pin will result in Sr equal to 
the identity matrix, so Sr can be omitted. 



Transform 



New 



(/ + e,<(N/^-l))t/,, where U,=\I-2^ 
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[0085] 
[0086] 
[0087] 



Theoretical Basis for Square Root Factorization 
Theory for Generalized QR Factorization 

The U defined in the generalized QR factorization has the property 



that 



( 










I 




S 


I 














V 











u7sis.. = 



[0088] U also has the property that 



1-2 



v-v S * 



2-v'"5-x 



'x!S x,+S,(t' -l a xje, 



/ ^2(xJSXi-(T-xJei) 



2-x;S-Xi -2-(T-x;ei 

where xi was used in defining v. 

[0089] Thus, the transform zeros out all the elements of a column matrix but 
one, which is set to c. The use of this later property to transform a matrix into triangular 
form is well known. When this is done there is a sequence of transforms taking Pin to P: 
P = [(/ + e,ej {J^ - 1)). U, J- [(/ + eje] {j^ -l}u j\- [(/ + e,el - 1)> U, \pin 

[0090] Due to the structure defined for U and e in the sequence, this can be 
rearranged into the form 

P=S'^-U-Pin where 

Sl^ ^ t + e,ej (J^ - 1)). (/ + ejc] (j^ - 1)) (/ + e,el (j^ - 1)) 



u=u,-u, 



[0091] Sr and Sw are related S;"*S"' = 5^," . Replacing Sw with Sr yields 
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P^SgP = Pin^U^S • U • Pin = Pin^S • Pin 
[0092] where the last equation is the result of multiplying each side of the 

previous on its left by its respective transpose. The last equation shows this is the desired 

transform. 

[0093] Theory for Factoring H 

[0094] The theory of the factorization of H is based on identities formed by 
multiplying each side of Eq. 1 by its respective transpose for each time point to yield 

M'' -M+G^ -G-A^ W M A = C''Q^''^Q"''C 
F^ ■G-B' M A = D'Q^'^Q"^C 

[0095] These identities are known in control theory textbooks. Their use to 
factor H and LSM is novel. The factorization of H can be verified by symbolically 
expanding the expression W 'L^SL ' using the definitions for L, S, and W in Eqs. 2 
and 3., and then using these identities. 

[0096] Theory for Factoring the LSME for IP Method 
The factorization in Eq. 5 implies 

"/ 



[0097] 



p^5 P=[r^r'' 



u 



L 



L 



= }^T-'T-^J + VSL 



Rearrange this to get the needed relation 
W^P^SPW-^ -J''T-'T-^J = W*VSL-W-^ =H 

[0098] If we substitute Eq 8 into Eq 7 and use the fact that S ' = S, we get 





'S 


P 0' 


mz 




'w-f 


0 


-I 


_J-Tj J 


m 




K 
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[0099] Scale the first block of rows by W' and substitute in for mz to yield 





'S 


'p-w^ o' 


z 




'f 
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-I 


_J-Tj J 


m 




K 



[00100] Symbolic perfonning the matrix multiplications 



J -r'"rj[mj \k_ 



z 




7' 


m 




K 



H r 



[00101] This shows the Square Root Algorithm for IP solves the LSME, 
presuming the algorithm for performing generalized QR factorization of Eq. 5 works. 
[00102] Theory for Factoring the LSME for AS Method 
[00103] The factorization in Eq. 9 implies 



=> P ' SgP = J ■ L-'S • L-y ' , where S„ = Sl'^S ■ S^^ 
[00104] If we substitute Eq 12 into Eq 1 1 and use the fact that S ' = S, we get 



P"^ JU'S 
0 



■-S. OT P 01 

_ 0 SJ[5L-^J'^ lJL, 



m 
mz 





K 




W-f 



[00105] Scale the second block of rows by W and substitute in for mz to yield 
P"^ JU^S 

0 w-'iJ 

[00106] Symbolic performing the matrix multiplications 
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m 




K 
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J* 

- J \\m 
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[00107] This shows the Square Root Algorithm for AS solves the LSME, 
presuming the algorithm for performing generalized QR factorization of Eq. 9 works. 
[00108] Computational Efficiency of the Square Root Method 
[00109] . The advantages of this factorization are that 

[00110] The LSME matrix is factored into the product of upper triangular, 
lower triangular matrix and diagonal matrices. 

[00111] The LSME can then be solved with a block back substitution pass, 
with the block upper triangular matrix followed by a forward substitution pass with the 
block lower triangular matrix. 

[00112] LSM and H do not need to be determined. 

[00113] The computer wordlength needed for L, M, and P, is only half that 
needed for H and LSM for a given level of numerical accuracy. The ability to use single 
precision computer arithmetic, instead of double precision, approximately doubles the 
computation speed. 

[00114] The factorization of H needs only be done once per controller update. 
Each iteration of the QP solver reuses this factorization. 

[00115] Completion of the factorization each QP iteration exploits the banded 
structure. Fully for use in IP methods; partially in AS methods. 

[00116] Further numerical efficiencies can be achieved with special QR and 
triangular matrix solvers that exploit the remaining structure of L and P. 

[00117] This algorithm for factoring H has the following computational 
advantages: 

[00118] The M and F matrices are lower triangular, which greatly reduces the 
computation in matrix operations with them, especially when N gets large. 
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[00119] The L and S matrices are sparse, banded matrices, which also reduces 
computations 

[00120] Many of the "square root" methods of control theory can be applied. 
Some of these exploit the structure within matrices A, B, C, D, Q, and R to further reduce 
computation. 

[00121] The amount of computation required to factor H grows only linearly 
with N, and 

[00122] The repeating structures of H and J are exploited. 

[00123] Figure 1 is a schematic of a model predictive control system 10 
according to the present invention for performing the algorithm described above. The 
system 10 controls a plant 12 with actuators 14 based upon feedback from sensors 16 and 
based upon commands. The system 10 includes a desired trajectory generator 20 for 
generating the desired trajectory based upon the commands. The system 10 includes a 
quadratic programming module 21 including a module 22 for determining the various 
matrices and vectors of the KKT conditions in the manner described above. The KKT 
conditions are also functions of the current state of the system 10 based upon feedback 
from sensors 16. The quadratic programming module 21 further includes a module 24 
for determining Hr as described above. A Modified Quadratic Programming Solver 26 
follows the algorithm described above to generate actuator conmiands to the actuators 14 
to bring about the desired state of the plant 12 and system 10. 

[00124] Figure 2 illustrates the standard way quadratic progranraiing (QP) 
solvers are embedded in a Model Predictive Controller (MPC) and how the novel square 
root method is embedded into standard QP solvers. The algorithm of Figure 2 is 
executed every controller update. Steps 50 and 52 of the figure translate model 
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predictive control into a quadratic programming problem and extract the optimal 
trajectories for the next N controller updates, including the current actuator conmiand, 
from the QP solution, respectively. All other steps are associated with the QP solver. 
General QP solvers consist of a structured, iterative search for the optimal solution that 
meets constraints. This consists of analyzing the solution from the prior iterate and setting 
up a version of the LSME (although not generally in the structured form presented here) 
in step 56, solving the LSME, and analyzing the results to determine if the solution has 
been achieved to the required accuracy in step 62. As shown in Figure 2, the square root 
method adds a new logic to precompute Hr in step 54, the square root of H, before the 
iterative search is begun. It also replaces the solution of the LSME with a novel 
computationally more efficient algorithm in step 60, including forming LSMroot in step 
58. 

[00125] The square root method also changes the interface between the top 
block of Figure 2 and the QP solver. In step 50, the commands for the system under 
control are received, a desired response to those commands is formulated, sensors signals 
indicating the current state of the system are received, and a QP problem as per the KKT 
conditions is formulated. The equations of the KKT condition are generally gathered into 
three matrix equations: one associated with the performance index, one associated with 
system dynamics model, and one associated with the limit conditions. The prior solution 
for the optimal trajectories may be used in forming the various matrices in the KKT 
conditions. In the standard method, the QP then forms a version of the LSME, with 
different structure, each iteration. In contrast, the square root method also modifies that 
portion of the QP solver, in step 56, in that the various variables, equations and matrices 
of the KKT conditions are gathered into the LSME in a specific order and structure. 
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[00126] In accordance with the provisions of the patent statutes and 
jurisprudence, exemplary configurations described above are considered to represent a 
preferred embodiment of the invention. However, it should be noted that the invention 
can be practiced otherwise than as specifically illustrated and described without departing 
from its spirit or scope. 
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